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ABSTRACT 

Aims. Two new families of self-consistent axisymmetric truncated equilibrium models for the description of quasi-relaxed rotating 
stellar systems are presented. The first extends the well-known spherical King models to the case of solid-body rotation. The second 
is characterized by differential rotation, designed to be rigid in the central regions and to vanish in the outer parts, where the imposed 
energy truncation becomes effective. 

Methods. The models are constructed by solving the relevant nonlinear Poisson equation for the self-consistent mean-field potential. 
For rigidly rotating configurations, the solutions are obtained by an asymptotic expansion based on the rotation strength parameter, 
following a procedure developed earlier by us for the case of tidally generated triaxial models. The differentially rotating models are 
constructed by means of a spectral iterative approach, with a numerical scheme based on a Legendre series expansion of the density 
and the potential. 

Results. The two classes of models exhibit complementary properties. The rigidly rotating configurations are flattened toward the 
equatorial plane, with deviations from spherical symmetry that increase with the distance from the center. For models of the second 
family, the deviations from spherical symmetry are strongest in the central region, whereas the outer parts tend to be quasi-spherical. 
The relevant parameter spaces are thoroughly explored and the corresponding intrinsic and projected structural properties are de- 
scribed. Special attention is given to the effect of different options for the truncation of the distribution function in phase space. 
Conclusions. Models in the moderate rotation regime are best suited to applications to globular clusters. For general interest in stellar 
dynamics, at high values of the rotation strength the differentially rotating models tend to exhibit a toroidal core embedded in an oth- 
erwise quasi-spherical configuration. Physically simple analytical models of the kind presented here provide insights into dynamical 
mechanisms and may be a useful basis for more realistic investigations carried out with the help of N-body simulations. 

Key words. Methods: analytical - globular clusters: general 



1. Introduction 

A large class of stellar systems is expected to be in a quasi- 
equilibrium state characterized by a distribution function not far 
from a Maxwellian. Mechanisms that are thought to operate in 
this direction range from standard two-body gravitational scat- 
tering (Chandrasekhar 1943 ) to violent relaxation, instabilities, 
and phase mixing for collisionless systems (Lynden-Bell 1967). 
In particular, for relatively small stellar systems, such as glob- 
ular clusters, two-star collisional relaxation is thought to be re- 
sponsible for such quasi-relaxed condition, while for large stel- 
lar systems, such as elliptical galaxies, incomplete violent re- 
laxation has likely acted during the formation process so as to 
bring the system to a state in which significant (radially biased) 
anisotropy is present only for weakly bound stars. To be sure, 
full relaxation would lead to a singular state, characterized by 
infinite mass; the so-called isothermal sphere solution may pro- 
vide an approximate representation only of the inner parts of real 
stellar systems. 

In line with the above considerations, several physically mo- 
tivated self-consistent finite-mass models (thus characterized by 
significant, but not arbitrary, deviations from a pure Maxwellian) 
have been constructed and studied, with interesting astronom- 
ical applications. For the description of partially relaxed sys- 
tems, self-consistent anisotropic models have been introduced 
to describe the products of incomplete violent relaxation (see 
Bertin & Stiavelli 1993 and Trenti et al. [2)05; and references 
therein) . For the description of quasi-relaxed globular clusters, 



the spherical isotropic King (1966) models, which incorporate 
in a heuristic way the concept of tidal truncation, provide a stan- 
dard reference framework (see Zocchi et al. 120121 for a photo- 
metric and kinematic study of a sample of Galactic globular 
clusters under different relaxation conditions). Recently, King 
models have been extended to the full triaxial case, to include 
explicitly the three-dimensional role of an external tidal field in 
the simple case of a circular orbit of the cluster within the host 
galaxy (Bertin & Varri 120081 Varri & Bertin 120091 in the fol- 
lowing denoted as Paper I and II, respectively; see also Heggie 
& Ramamani 1995); these triaxial models can thus be seen as a 
collisionless analogues of the purely tidal Roche ellipsoids. 

For astronomical applications, one of the main drivers (or 
empirical clues) in the above modeling investigations is the is- 
sue of the origin of the observed geometry of stellar systems, 
being it known that pressure anisotropy, tides, or rotation can be 
responsible, separately, for deviations from spherical symmetry. 
In the above-mentioned studies, very little attention was actu- 
ally placed on the role of rotation. For ellipticals, most of the 
attention that led to the development of stellar dynamical mod- 
els, after the first kinematical measurements became available 
in the mid-70s, was taken by the study of the curious behav- 
ior of pressure-supported systems in the presence of anisotropic 
orbits (see also Schwarzschild [19791 [19521 de Zeeuw [l985l In 
contrast, very little effort has been made in modeling rotation- 
dominated ellipticals, even though the entire low-mass end of the 
distribution of elliptical galaxies might be consistent with a pic- 
ture of rotation-induced flattening (e.g., see Davies et al. 1983 
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and more recently Emsellem et al. 1201 II) ; similar comments ap- 
ply to bulges. 

For globular clusters, given the fact that they only exhibit 
modest amounts of flattening and given the success of the spher- 
ical King models, little work has been carried out in the direction 
of stationary self-consistent rotating models (with some notable 
exceptions, that is Wolley & Dickens 1962, Lynden-Bell 1962, 
Kormendy & Anand 119711 Lupton & Gunn 119871 and Lagoute 
& Longaretti 1996 ). Therefore, as far as rotation-dominated sys- 
tems are concerned, much of the currently available modeling 
tools go back to the pioneering work of Prendergast & Tomer 
d 19701 ). Wilson (119751 ). and Toomre d 19821), intended to d escrib e 
ellipticals, and of Jarvis & Freeman (119851 ) and Rowley ( 1988), 
devoted to bulges. In general, we may say that only very few 
rotating models with explicit distribution function are presently 
known (for a recent example, see Monari et al. in preparation). 
In this context, one should also mention the interesting work 
by Vandervoort (1980) on the collisionless analogues of the 
Maclaurin and Jacobi ellipsoids. 

On the empirical side, a deeper study of quasi-relaxed ro- 
tating stellar systems is actually encouraged by the widespread 
conviction that the geometry of the inner parts of globular clus- 
ters, where some flattening is noted, should not be the result of 
tides, but rather of rotation (King 1 1961b . In other words, it is 
frequently believed that tides and pressure anisotropy (and dust 
obscuration), even though playing some role in individual cases, 
should not be considered as the primary explanation of the ob- 
served flattening of Galactic globular clusters. Such conclusion 
is suggested by the White & Shawl ( 1987 ) database of elliptici- 
ties. However, for this database we note that: (i) the cluster flat- 
tening values do not all refer to a standard isophote, such as the 
cluster half-light radius (as also noted by van den Bergh 2008), 
(ii) the data mostly refer to the inner regions. Such limitations 
are crucial because there is observational evidence that the el- 
lipticity of a cluster depends on radius (see Geyer et al. I1983I ). 
In turn, recent studies by Chen & Chen (120101 ). in which these 
restrictions are addressed more carefully, suggest that Galactic 
globular clusters are less round than previously reported, espe- 
cially those in the region of the Galactic Bulge; at variance with 
previous analyses, it seems that their major axes preferentially 
point toward the Galactic Center, as naturally expected if their 
shape is of tidal origin. 

The connection between flattening and internal rotation has 
been discussed in detail by means of nonspherical dynamical 
models in just a handful of cases, in particular for to Cen (for 
an oblate rotator nonparametric model, see Merritt 1997; for 
an orbit-based analysis, see van de Ven 120061 for an applica- 
tion of the Wilson [19731 mo dels, see Sollima et al. |2059l), 47 
Tuc (Meylan & Mavor lT986| M15 (van den Bosch et al. 2006 ), 
and M13 (Lupton et al. 1987 ). In addition, specifically designed 
2D Fokker-Planck models (Fiestas et al. 2006) have been ap- 
plied to the study of M5, NGC 2808, and NGC 5286. We re- 
call that the detection of internal rotation in globular clusters 
is a challenging task, because the typical value of the ratio of 
mean velocity to velocity dispersion is only of a few tenths, for 
example V/cr ~ 0.26,0.32 for 47 Tuc and a> Cen, respectively 
(Meylan and Mavor [T986t for a summary of the results for sev- 
eral Galactic objects, see also Table 7.2 in Meylan & Heggie 
119971 ). However, great progress made in the acquisition of pho- 
tometric and kinematical information, and in particular of the 
proper motion of thousands of stars (for to Cen, see van Leeuwen 
et al. 120001 Anderson & van der Marel 2010; for 47 Tucanae, see 
Anderson & King 2003 and McLaughlin et al. 120061 ). makes this 
goal within reach (see Lane et al. 2009, Lane et al. l2010l for new 



kinematical measurements, in which rotation, when present, is 
clearly identified). 

On the theoretical side, two general questions provide fur- 
ther motivation to study quasi-relaxed rotating stellar systems. 
On the one hand, many papers have studied the role of rotation 
in the general context of the dynamical evolution of globular 
clusters, but a solid interpretation is still missing. Early inves- 
tigations (Agekian 1958, Shapiro & Marchant 1976) suggested 
that initially rotating systems should experience a loss of angular 
momentum induced by evaporation, that is, angular momentum 
would be removed by stars escaping from the cluster. Because 
of the small number of particles, N-body simulations were ini- 
tially (Aarseth 119691 Wielen 119741 and Akiyama & Sugimoto 
1989) unable to clearly describe the complex interplay between 
relaxation and rotation. Later investigations, primarily based on 
a Fokker-Planck approach (Goodman 1983 Einsel & Spurzem 
[19991 Kim et al. 125521 and Fiestas et al. 125561) have clarified this 
point, not only by testing the proposed mechanism of angular 
momentum removal by escaping stars, but also by showing that 
rotation accelerates the entire dynamical evolution of the sys- 
tem. More recent N-body simulations (Boily 2000 Ernst et al. 
120071 and Kim et al. 120081 ) confirm these conclusions and show 
that, when a three-dimensional tidal field is included, such ac- 
celeration is enhanced even further. The mechanism of angular 
momentum removal is generally considered to be the reason why 
Galactic globular clusters are much rounder than the (younger) 
clusters in the Magellanic Clouds, for which an age-ellipticity re- 
lation has been noted (Frenk & Fall 1982), but other mechanisms 
might operate to produce the observed correlations (Meylan & 
Heggie [19971 van den Bergh 120081 . 

On the other hand, the role of angular momentum during the 
initial stages of cluster formation should be better clarified. In 
the context of dissipationless collapse, relatively few investiga- 
tions have considered the role of angular momentum in numeri- 
cal experiments of violent relaxation (e.g., the pioneering studies 
by Gott |19731 see also Aguilar & Merritt |T990l ). Interestingly, the 
final equilibrium configurations resulting from such collisionless 
collapse show a central region with solid body rotation, while the 
external parts are characterized by differential rotation. 

In conclusion, mastering the internal structure of spheroidal 
and triaxial stellar systems through a full spectrum of models, in- 
cluding rotation, is a prerequisite for studies of many empirical 
and theoretical issues. In addition, it is required for the inter- 
pretation of the relevant scaling laws (such as the Fundamental 
Plane, which appears to extend from the brightest, pressure sup- 
ported ellipticals down to the low-luminosity end of the distribu- 
tion of early-type galaxies, and possibly further down to the do- 
main of globular clusters) and for investigations aimed at iden- 
tifying the presence of invisible matter (in the form of central 
massive black holes or diffuse dark matter halos) from stellar 
dynamical measurements. It would thus be desirable to construct 
rotating models to be tested on low luminosity ellipticals, bulges, 
and globular clusters, especially now that important progress has 
been made in the collection and analysis of kinematical data. 
Presumably, many of these stellar systems are quasi-relaxed. In 
which directions should we explore deviations from the strictly 
relaxed case? 

In the present paper, we consider two families of axisym- 
metric rotating models: the first one is characterized by the pres- 
ence of solid-body rotation and isotropy in velocity space (see 
Appendix B in Paper I for a brief introduction). Indeed, full re- 
laxation in the presence of nonvanishing total angular momen- 
tum suggests the establishment of solid-body rotation through 
the dependence of the relevant distribution function / = f(H) 
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on the Jacobi integral H — E — a>J z (see Landau & Lifchitz 
119671 p. 125). But, for applications to real stellar systems, one 
may take advantage of the fact that the collisional relaxation time 
may be large in the outer regions, so that in the outer parts the 
constraint of solid-body rotation might be released. In particular, 
for globular clusters we may argue that the outer parts fall into 
a tide-dominated regime, for which evaporation tends to erase 
systematic rotation even if initially present, as confirmed by the 
above-mentioned study based on the Fokker-Planck method. For 
the truncation, we may then consider a heuristic prescription to 
simulate the effects of tides, much like for the spherical King 
models. 



In view of possible applications to globular clusters, we 
thus consider a second class of axisymmetric rotating models 
based on a distribution function dependent only on the energy 
and on the z-component of the angular momentum / = f(I) 
where / = I(E, J z ), with the property that I ~ E for stars with 
relatively high z-component of the angular momentum, while 
I ~ H - E - <±>J Z for relatively low values of J z . Such models 
are indeed defined in order to have differential rotation, designed 
to be rigid in the center and to vanish in the outer parts, where 
the energy truncation becomes effective. As far as the velocity 
dispersion is concerned, this family may show a variety of pro- 
files (depending on the values of the relevant free parameters), 
all of them characterized by the presence of isotropy in the cen- 
tral region. We thus add two classes of self-consistent models 
to the relatively short list of rotating stellar dynamical models 
currently available. 



One aspect that plays an important role in defining a phys- 
ically motivated distribution function, which often goes unno- 
ticed (but see Hunter [19"771 Davoust [T9T7l and Rowley [19881) . is 
the choice of the truncation prescription in phase space. The ad- 
vantages and the limitations of alternative options available for 
the second family of models will be discussed in detail. In this 
context, we will also address the issue of whether these differen- 
tially rotating models fall within the class of systems for which 
rotation is constant on cylinders. 



The paper is structured as follows. The properties of the fam- 
ily of rigidly rotating models, constructed on the basis of general 
statistical mechanics considerations, are illustrated in Sect. [2] 
The family of differentially rotating models, designed for the 
application to rotating globular clusters, is introduced in Sect. [3] 
where we briefly describe the method used for the solution of 
the self-consistent problem and discuss the relevant parameter 
space. Section|4]is devoted to a study of the intrinsic properties 
and Sect. [5] to the projected observables derived from differen- 
tially rotating models. After illustrating in Sect. [6] the effect of 
different truncation prescriptions in phase space, we summarize 
the results and present our conclusions in Sect. [7] The appen- 
dices are devoted, respectively, to a discussion of the nonrotat- 
ing limit of our families of models, to the details of the alterna- 
tive truncation option for the second family, and to a description 
of the code used for the construction of the differentially rotat- 
ing configurations. A study of the dynamical stability and of the 
long-term evolution of the families of models introduced here 
will be addressed by means of an extensive survey of specifically 
designed N-body simulations and will be presented in following 
papers (Varri et al. 1201 II Varri et al., in preparation). 




Fig. 1. The boundary surface, defined implicitly by \ji(r) - 0, 
of a critical second-order rigidly rotating model with *F = 2. 
The configuration is axisymmetric; the points on the equatorial 
plane are saddle points (see Sect. 2.2 for details). Rotation takes 
place around the z-axis. The spatial coordinates are expressed in 
appropriate dimensionless units (see Eq. ©). 



2. Rigidly rotating models 

2. 1 . The distribution function 

The construction of rigidly rotating configurations characterized 
by nonuniform density is a classical problem in the theory of 
rotating stars, starting with Milne dl923b and Chandrasekhar 
( I1933I ), but it basically remained limited to the study of a fluid 
with polytropic equation of state, for which the solution of 
the relevant Poisson equation can be obtained by means of a 
semi-analytical approach (for a comprehensive description, see 
Chaps. 5 and 10 in Tassoul [19781 for an enlightening presenta- 
tion of the general problem of rotating compressible masses, see 
Chap. 9 in Jeans [T9291 . The reader is referred to Paper I for a 
discussion of the application of some of the mathematical meth- 
ods developed in that context to the construction of nonspherical 
truncated self-consistent stellar dynamical models. The devia- 
tions from spherical symmetry studied in Paper I are induced 
by the presence of a stationary perturbation characterized by a 
quadrupolar structure, that is, either an external tidal field or in- 
ternal solid-body rotation. In particular, in Appendix B of Paper 
I we briefly outlined the extension of the King ( 1966) models to 
the case of internal solid-body rotation, of which, after a short 
summary of the relevant definitions, we now provide a full de- 
scription in terms of the relevant intrinsic and projected proper- 
ties. 

The extension of the family of King models is performed by 
considering the distribution function: 

f K (H) = Ae- aHo ^- a(H - Ho) - l] (1) 

for H < Hq and f' K {H) = otherwise, where 

H = E-coJ z (2) 

denotes the Jacobi integral, with a> the angular velocity of the 
rigid rotation (hence the superscript r), assumed to take place 
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around the z-axis. The quantities E and J, are the specific one- 
star energy and z-component of the angular momentum, Hq rep- 
resents a cut-off constant of the Jacobi integral, while A and a 
are positive constants. In the corotating frame of reference the 
Jacobi integral can be written as the sum of the kinetic energy, 
the centrifugal potential Q> cen {x,y) = —{x 2 + y 2 )u> 2 /2 (with the 
equatorial plane (x, y) perpendicular to the rotation axis given by 
the z-axis), and the cluster mean-field potential <&c, to be de- 
termined self-consistently. Therefore, the dimensionless escape 
energy can be expressed as: 



i/r(r) = a{H - [<l>c(r) + <S> cm {x,y)]} 



(3) 



and the boundary of the cluster, implicitly defined by i[r(r) = 0, 
is an equipotential surface for the total potential i>c + $>cen- Note 
that, by construction, in the limit of vanishing internal rotation, 
this family of models reduces to the family of spherical King 
( 1966) models (see AppendixlAlfor a summary of the main prop- 
erties of the family in the nonrotating limit). 

The construction of the models requires the integration of 
the associated nonlinear Poisson equation, which, after scaling 
the spatial coordinates with respect to the scale length 



ro = 



4nGpoa 



1/2 



can be written as 

p K (>t> m ) 



vy in,) = -9 

where 

= ^ 



-2* 



(4) 



(5) 



(6) 



is the dimensionless parameter that characterizes the rotation 
strength and 



(7) 



is the depth of the potential well at the center. The dimensionless 
density profile is given by 



where 

8tt2 1/2 A 



A = 



3a 3 / 2 



-aH 



(8) 



(9) 



and y denotes the incomplete gamma function. Therefore, the 
central density is given by po = Apx CP). Outside the cluster, for 
negative values of if/, we should refer to the Laplace equation 



VV ( "° = 18*. 



(10) 



The relevant boundary conditions are given by the require- 
ment of regularity of the solution at the origin and by the con- 
dition that + a$> C en — * oHo at large radii. The Poisson 
(internal) and Laplace (external) domains are thus separated by 
the boundary surface, which is unknown a priori. Therefore, we 
have to solve an elliptic partial differential equation in a free 
boundary problem. In particular, here we illustrate the proper- 
ties of the solutions of the Poisson-Laplace equation obtained 
by using a perturbation method, which also requires an expan- 
sion of the solution in Legendre series. To obtain a uniformly 




Fig. 2. Parameter space for second-order rigidly rotating mod- 
els. The solid line represents the critical values of the rotation 
strength parameter % cr and the grey region identifies the values 
QV,X) f° r which the resulting models are bounded by a closed 
constant-t/' surface (subcritical models). 

valid solution over the entire space, an asymptotic matching is 
performed between the internal and the external solution, using 
the Van Dyke principle (see Van Dyke 1 1 975b . This method of 
solution is basically the same as proposed by Smith d 19751 1 for 
the construction of rotating configurations with polytropic index 
n = 3/2. We have calculated the complete solution up to second- 
order in the rotation strength parameter x- The final solution is 
expressed in spherical coordinates f = (f, 8, </>) and the resulting 
configurations are characterized by axisymmetry (i.e., the den- 
sity distribution and the potential do not depend on the azimuthal 
angle </>). For the details of the method for the construction of the 
solution the reader is referred to Paper I and to Varri (120121 1. in 
which the complete calculation is provided. 

2.2. The parameter space 

The resulting models are characterized by two dimensional 
scales (e.g., the total mass and the core radius) and two di- 
mensionless parameters. As in the spherical King models, the 
first parameter measures the concentration of the configuration; 
we thus consider the quantitjQ (see Eq. ©) or equivalently 
c = log(r fr /ro), where r tr is the truncation radius of the spherical 
King model associated with a given value of T*. The second di- 
mensionless parameter* (see Eq. (O) characterizes the rotation 
strength measured in terms of the frequency associated with the 
central density of the cluster. 

For every value of the dimensionless central concentration 
m there exists a maximum value of the rotation strength param- 
eter, corresponding to a critical model, for which the boundary 
is given by the critical constant-i/r surface. The boundary surface 
of a representative critical model (with *P = 2) is depicted in 
Fig. [1] the surface is such that all the points on the equatorial 
plane are saddle points, where the centrifugal force balances the 
self-gravity. We refer to their distance from the origin as rg, the 



1 In the literature, the parameter *P is often denoted by Wo; here we 
prefer to keep the same notation used in Paper I, Paper II, and in other 
articles. 
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break-off radius. In the constant-^ family of surfaces associated 
with a given value of "F, the critical surface thus separates the 
open from the closed surfaces. Consistent with the assumption 
of stationarity, only configurations bounded by closed surfaces 
are considered here. This unique geometrical characterization 
suggests that the effect of the rotation may be expressed also 
in terms of an extension parameter 

6=^, (ID 
rB 

which provides an indirect measure of the deviations from 
sphericity of a configuration, by considering the ratio between 
the truncation radius of the corresponding spherical King model 
and the break-off radius of the associated critical surface. 
Therefore, a given model may be labelled by the pair Q¥,x) or 
equivalently by the pair Q¥, 6). For given T*, there is thus a maxi- 
mum value of the allowed rotation, which we may express a&Xcr 
or 5 cr . A model with 6 < 6 cr may be called subcritical. 

For each value of W, the critical value of the rotation param- 
eter can be found by numerically solving the systerrQ 

(d ft f f (r = h,e = 7r/2;xcr)=0 

\^(f=f B ,e = it/2;Xcr)=0, ' " ] 

where the unknowns are fs and Xcr- In terms of the extension 
parameter, for a given *P, the critical condition occurs when 6 = 
6 cr = r, r /rB ~ 2/3. This value is obtained by inserting in Eq. (TTZb 
the zeroth-order expression for the cluster potential, as discussed 
in detail in Sect. 2 of Paper II for the tidal problem (see p. 251- 
255 in Jeans 1929 for an equivalent discussion referred to the 
purely rotating Roche model, that is a rotating configuration in 
which a small region with infinite density is surrounded by an 
"atmosphere" of negligible mass). 

The parameter space for the second-order models is pre- 
sented in Fig.|2](which corresponds to Fig. 1 of Paper II describ- 
ing the tidal models). Two rotation regimes exist, namely the 
regime of low-deformation (6 <K S cr , bottom left corner), where 
internal rotation does not affect significantly the morphology of 
the configuration, which remains very close to spherical symme- 
try, and that of high-deformation (6 ~ S cr , close to the solid line), 
where the model is highly affected by the nearly critical rotation 
velocity, especially in the outer parts. Note that the actual regime 
depends on the combined effect of rotation strength and of con- 
centration. In other words, the models described here belong to 
the class of rotating configurations characterized by equatorial 
break-off ("region of equatorial break-off", see Fig. 44, p. 267 in 
Jeans 1 19291 1, for which the limiting case is given by the purely 
rotating Roche model. 

A global kinematical characterization, complementary to the 
information provided by the rotation strength parameter, is of- 
fered by the parameter t — K or d/\W\, defined as the ratio between 
ordered kinetic energy and gravitational energy. Figure [3] illus- 
trates the relation between the two parameters for models with 

2 For brevity, here we omit the explicit structure of the system. The 
reader is referred to Eqs. (8)-(14) in Paper II, with respect to which sev- 
eral differences occur, because in the axisymmetric rotation problem the 
coefficients of the asymptotic series are best expanded in (normalized) 
Legendre polynomials (rather than spherical harmonics). In particular, 
with respect to Paper II: (i) the coefficients with m j= must be dropped; 
(ii) the coefficients with / 1 must be multiplied by the factor (n/4) </2 , 
because of the different normalization used in the two systems of or- 
thonormal functions; (iii) the expressions are evaluated at r B instead of 
f T . As in Paper I, for the definition of the Legendre polynomials we 
refer to Eqs. (22.3.8) and (22.2.10) of Abramowitz & Stegun ( fT965V 




0.0 0.2 0.4 0.6 0.8 1.0 

X/Xcr 



Fig. 3. Values of the ratio between ordered kinetic energy and 
gravitational energy t = K orc {/\W\ for selected rigidly rotating 
models characterized by *F = 1,3,5,7 (solid lines, from top to 
bottom) and^ in the range [0,^ cr ]. For comparison, the dashed 
line indicates the values of t for the sequence of Maclaurin oblate 
spheroids in the limit of small eccentricity. 



selected values of ¥ and increasing values of x, up to the criti- 
cal configuration characterized by Xcr- The parameter t increases 
linearly for increasing values of x (with a slope dependent on 
the concentration parameter *F). A similar linear relation is ob- 
served for small values of eccentricity (e « 1) in the sequence 
of Maclaurin oblate spheroids t(e) ~ x( e ) ~ 2e 2 /l5 (for the 
definitions of the two parameters in the context of Maclaurin 
spheroids, see Eqs. (10.20) and (10.24) in Bertin l2000l ; in Fig. 
[3] the relevant linear relation is normalized with respect to the 
maximum value of the rotation strength parameter attained in the 
sequence of Maclaurin spheroids Xmax = 0.11233 (see Eq. (10), 
p. 80 in ChandrasekMr [l969l . 

2.3. Intrinsic properties 

The geometry of the models, reflecting the properties of the cen- 
trifugal potential, is characterized by symmetry around the z- 
axis and reflection symmetry with respect to the equatorial plane 
(x,y). As expected, compared to the corresponding spherical 
King models, the rotating models stretch out on the equatorial 
plane and are slightly flattened along the direction of the rotation 
axis (see Fig.|4]for the density profiles of selected critical second 
order models evaluated on the equatorial plane and along the z- 
axis). In general, configurations in the low-deformation regime 
(6 <K 6 cr ), regardless of the value of concentration *P, are almost 
indistinguishable from the corresponding spherical King mod- 
els. 

Models in the intermediate and high-deformation regime 
(6 ~ 6 cr ) show modest deviations from spherical symmetry 
in the central regions, while they are significantly flattened in 
the outer parts. The intrinsic eccentricity profile, defined as 
e = [1 — (b/a) 2 ] 1 / 2 , where a and b are semi-major and semi- 
minor axes of the isodensity surfaces, is a monotonically increas- 
ing function of the semi-major axis (see Fig.|5]for the computed 
eccentricity profiles of selected critical second-order models). 
We recall that the geometry of the boundary surface of a crit- 



5 



A. L. Varri and G. Bertin: Self-consistent models of quasi-relaxed rotating stellar systems 




Fig. 4. Intrinsic density profiles (normalized to the central value) 
evaluated on the equatorial plane (solid lines) and along the z- 
axis (dashed lines) for critical second-order rigidly rotating mod- 
els with >P = 1,2, ...,10 (from left to right). 



Fig. 5. Eccentricity profiles of the isodensity surfaces for se- 
lected critical second-order rigidly rotating models, with *F = 
1,3,5,7 (from top to bottom); dotted horizontal lines show the 
central eccentricities values estimated analytically (see Eq. (fT4l>). 



ical model depends only slightly on the value of the concen- 
tration parameter. In particular, in the critical case, the break- 
off radius rg represents the distance from the center of the out- 
ermost points of the boundary surface on the equatorial plane 
and the truncation radius r, r is approximately the distance of 
the last point on the polar axis (i.e., the z-axis). Therefore, since 
&cr = ftrlrs ~ 2/3, the value of the termination points of the ec- 
centricity profiles of critical models is approximately the same 
(<? =s 0.75; see the termination points of the solid lines in Fig. [5]). 

In addition, by using the multipolar structure of the solution 
of the Poisson-Laplace equation obtained with the perturbation 
method, the asymptotic behavior of the eccentricity profiles in 
the central regions can also be evaluated analytically. Since the 
distribution function depends only on the Jacobi integral (i.e., 
the isolating energy integral in the corotating frame of refer- 
ence), the density and the velocity dispersion profiles are func- 
tions of only the escape energy (see Eqs. ([8]l and ( fTTl i, respec- 
tively). Therefore, there is a one-to-one correspondence between 
equipotential, isodensity, and isobaric surfaces and the eccentric- 
ity profiles can be calculated with reference to just one of these 
families of surfaces. In fact, if expanded to second order in the 
dimensionless radius, the escape energy in the internal region 
reduces to: 

f m {f, 9) = y ¥+^[-3+6x + 2A 2 U 2 (6) X + 

(B 2 + l)U 2 (6)x 2 ] r 2 + 0(r 4 ) , (13) 

where U 2 (0) denotes the normalized Legendre polynomial with 
I = 2 and A 2 , B 2 are appropriate (negative) coefficients, depend- 
ing on *P, which are determined by asymptotically matching 
the internal and external solution, in order to have continuity 
on the entire domain (see Eqs. (62) and (68) in Paper I, to be 
interpreted as indicated in Appendix B of Paper I). By setting 
if/(""\a, ji /2) = tf/(""\b, 0), we thus find that, in the innermost re- 
gion, the eccentricity tends to the following nonvanishing central 
value 

[6A 2X + 3(B 2 + l) X 2 ]V 2 

eo ■ . (14) 

[6 yfflS&X - 1) + \A 2X + 2(B 2 + l)x 2 ] l/z 



Therefore, the central value of the eccentricity is finite, of or- 
der 0(x l ^ 2 ), and strictly vanishes only in the limit of vanishing 
rotation strength. This result is nontrivial because the centrifu- 
gal potential (which induces the deviations from sphericity) is a 
homogeneous function of the spatial coordinates. Therefore, we 
might naively expect that, in their central regions, the models re- 
duce to a perfectly spherical shape (i.e., eo = 0), even for finite 
values of the rotation strength. This property has been noted also 
in the family of triaxial tidal models, in which the tidal potential 
plays the role of the centrifugal potential (see Sect. 3.1 of Paper 
II). 

Deviations from spherical symmetry can also be described, 
in a global way, by the quadrupole moment tensor, defined as 

Qij = I (3xiXj - r 2 6ij)p(r)d 2 r = 
Jv 

Ar 5 f {3x i x j -r 2 6 i j)p{r)d^ = ArlQ ij , (15) 
Jv 

where the integration is performed in the volume V of the entire 
configuration. It can be easily shown that in our coordinate sys- 
tem the tensor is diagonal and that Q xx = Q yy . The components 
of the tensor can be calculated explicitly 

m = QfJ = = 2*^8*0*0 (a* + ^) , (16) 

The quantities a 2 and b 2 are appropriate (positive) coeffi- 
cients, depending on T, resulting from the asymptotic match- 
ing of the internal and external solution of the Poisson-Laplace 
equation. The sign of the components are consistent with the 
above mentioned compression and stretching of the density dis- 
tribution. The expression in Eq. ( fTSI l is calculated from the 
second-order external solution; the first-order expression is re- 
covered by dropping the quadratic term in the parameter x- At 
variance with the tidal case, the ratio Q zz /Q xx is independent of 
the rotation parameter x- This analytical result has been com- 
pared to the ratio of the components of the quadrupole tensor 
determined by direct numerical integration (performed by means 
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Fig. 6. Intrinsic velocity dispersion profiles (normalized to the 
central value) evaluated on the equatorial plane (solid lines) and 
along the z-axis (dashed lines) for selected critical second-order 
rigidly rotating models with *F = 1,2, 10 (from left to right). 



of the algorithm VEGAS, see Press et al. 1 19921 ) and good agree- 
ment has been found. Q For the tidal case, the detailed calculation 
can be found in Sect. 3.3 and Appendix B of Paper II; the calcu- 
lation is easily adapted to the rotating case. 

By construction, the models are isotropic in velocity space, 
with the dimensionless scalar velocity dispersion given by 



2 y(7/2,<A) 
5 r(5/2, ^r) 



(17) 



As noted for the intrinsic density profiles, a compression along 
the vertical axis and a stretching along the equatorial plane occur 
also for the velocity dispersion profiles (the profiles of selected 
critical second-order models are shown in Fig. |6j>. 

The mean rotation velocity (which is subtracted away when 
the corotating frame of reference is considered) characterizing 
the models is defined as (v) = we.xc e,- = (v^)e^; in the adopted 
dimensionless units, the azimuthal component can be written as 



(v<p) = {vi>a 



1/2 



3x l/2 r sine. 



(18) 



As expected, the mean velocity is constant on cylinders (in our 
coordinate system, the cylindrical radius is defined by R — 
r sin 9). The relevant dimensionless angular velocity is linked to 
the rotation strength parameter by the following relation u> = 
3x 1 ^ 2 (the numerical factor 3 is due to the adopted scale length, 
see Eq. (|4|i). The rotation profiles of selected second-order crit- 
ical models are represented in Fig. [7] For all the models, as we 
approach the boundary of the configuration, the ratio (v^,)/&k 
quickly diverges since at the boundary the rotation velocity tends 
to a finite value while the velocity dispersion vanishes; this be- 
havior is observed in every direction (except for the z-axis, on 
which the rotation is absent by definition). 



3 Strictly speaking, the analytical expressions for the components of 
the quadrupole moment tensor refer only to the inner region, because 
the contribution from the boundary layer, where i]/ is of order 0(x), is 
neglected. 




Fig. 7. Mean rotation velocity profiles on the equatorial plane 
for the selected critical second-order rigidly rotating models with 
= 1,2, 10 (from left to right). The mean velocity increases 
linearly with radius because the rotation is rigid (see Eq. ( fT8l ; 
note the logarithmic scale of the horizontal axis). The values of 
the termination points of the curves depend on the value of Xcr 
(which decreases as *P increases, see Fig. |2]i and on the exten- 
sion of the models on the equatorial plane (which increases as *P 
increases, see Fig.|4|i. 



2.4. Projected properties 

For a comparison of the models with the observations (under 
the assumption of a constant mass-to-light ratio), we have then 
computed surface (projected) density profiles and isophotes. The 
projection has been performed along selected directions, identi- 
fied by the viewing angle (6, <p) corresponding to the zp axis of 
a new coordinate system related to the intrinsic system by the 
transformation Xp = Rx; the rotation matrix R = Ri(ff)RT,((f>) is 
expressed in terms of the viewing angles, by taking the xp axis 
as the line of nodes (i.e., the same projection rule we adopted 
for the triaxial tidal models, see Sect. 3.4 in Paper II). Since the 
rigidly rotating models are characterized by axisymmetry with 
respect to the z-axis and reflection symmetry with respect to 
the equatorial plane, it is sufficient to choose the viewing an- 
gles from the (Jc, z)-plane of the intrinsic coordinate system. In 
particular, we used the line of sights defined by 8/ = i(n/%) and 
= with i = 0, ..,4, and we calculated (by numerical inte- 
gration, using the Romberg's rule) the dimensionless projected 
density 



p(r P )dzp , 



(19) 



where 



■sp 



(x\ - xp - yp) 1/2 with x e the edge of the model 
along the x axis of the intrinsic coordinate system. The projec- 
tion plane (xp, yp) has been sampled on an equally-spaced square 
cartesian grid centered at the origin. 

The first five panels of Fig. [8] show the projected images of 
a critical second-order model with *F = 2; the first and the fifth 
panel correspond, respectively, to the least and to the most fa- 
vorable line of sight for the detection of the intrinsic flattening 
of the model (the z and x-axis of the intrinsic coordinate system, 
that is "face-on" and "edge-on" view). 
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Fig. 8. Projections along directions identified by <p — and 9 = /(tt/8) with i = 0, .., 4 (from left to right, top to bottom) of a critical 
second-order rigidly rotating model with *P = 2; the first and the fifth panel represent the projections along the z-axis ("face on") 
and the i-axis ("edge-on") of the intrinsic coordinate system, respectively. Solid lines mark the isophotes, corresponding to selected 
values of Z/Xo in the range [0.9, 10~ 7 ]. The last panel shows the ellipticity profiles, as functions of the semi-major axis of the 
projected image ap, referred to the lines of sight considered in the previous panels (from bottom to top). Dots represent the locations 
of the isophotes and the arrow marks the position of the half-light isophote (practically the same for every projection considered in 
the figure). 



The morphology of the isophotes of a given projected image 
can be described in terms of the ellipticity profile, defined as 
s = 1 - bpjcip where ap and bp are the principal semi-axes, 
as a function of the semi-major axis ap. As already noted for 
the (intrinsic) eccentricity profile, the deviation from circularity 
increases with the distance from the origin. In the inner region, 
the central value of the ellipticity is consistent with the central 
eccentricity <?o calculated in the previous subsection. The last 
panel of Fig.[8]illustrates the ellipticity profiles corresponding to 
the projections displayed in the previous panels. In addition, the 
isophotes of models in the high deformation regime (6 ~ 5 cr ), if 
projected along appropriate line of sights, show clear departures 
from a pure ellipse, that can be characterized as a "disky" overall 
trend (e.g., see Jedrzejewski 1987), particularly evident in the 
outer parts (see fourth and fifth panels in Fig. [8]). 



Because the family of rigidly rotating models is character- 
ized by simple kinematical properties (pressure isotropy and 
solid-body rotation), that have been already presented in detail 
with reference to three-dimensional configurations, for brevity, 
the derivation of the projected kinematical properties is omitted 
here; a full description can be found in Varri ( 120121) . 



3. Differentially rotating models 

3. 1 . Choice of the distribution function 

As indicated in the Introduction, theoretical and observational 
motivations have brought us to look for more realistic config- 
urations, characterized by differential rotation. Thus we focus 
our attention on axisymmetric systems, within the class of dis- 
tribution functions that depend only on the energy E and the 
z-component of the angular momentum J z , and we consider the 
integral 

I(E, J Z )=E , (20) 

1 + bj} c 

where to, b, and c > 1/2 are positive constants. The quantity 
I(E, J z ) reduces to the Jacobi integral for small values of the z- 
component of the angular momentum and tends to the single- 
star energy in the limit of high values of J z . Therefore, if we 
refer to a distribution function of the form / = /(/), we may 
argue that co is related to the angular velocity in the central re- 
gion of the system, characterized by approximately solid-body 
rotation, whereas the positive constants b, c will determine the 
shape of the radial profile of the rotation profile. In view of the 
arguments that have led to the truncation prescription that char- 
acterizes King model, we decided to introduce a truncation in 
phase-space based exclusively on the single-star energy with re- 
spect to a cut-off constant Eq. 
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Table 1. Summary of the properties of the families of models studied in the present paper. 



Family of 


Distribution 


Dimensionless 


Internal 


Anisotropy 


Isophote 


Paper 


models 


function 


parameters 


rotation 


profile 


shape 


sections 


fXP) 




%X 


solid-body 


(0, 0,0) 


disky 


2 




Ae - fl £„ [e - fl( ;-£ ) _ i + a (J_ £ )] 


*¥, X, b,c 


differential 


(0,> 0,-2) 


boxy 


3,4,5,6 




Aa-aEoo-aV-Eo) 


*¥,X, b,c 


differential 


(0,> 0,0) 


boxy 


3,6, App. B] 



Notes. The relevant integrals are defined as H = E - ljJ- and I = E - ojJ-J{\ + bJ: c ), with the corresponding cut-off constants given by H and 
Eq. The dimensionless parameters are defined as follows: *P = if/(0) respresents a measure of the concentration, x = u 2 /(AnGpo) a measure of the 
(central) rotation strength, and, for the family of differentially rotating models, b = brjfcr c and c > 1/2 determine the shape of the rotation profile. 
The pressure anisotropy profiles are characterized in terms of the values of the anisotropy parameter a = 1 — o~^/o~ 2 rr in the central, intermediate, 
and outer regions of a model; values of a greater than, lower than, and equal to zero indicate radially-biased, tangentially-biased anisotropy, and 
isotropy in velocity space, respectively. For each family of models, the first and third values of a are calculated analytically as the limiting values 
for small and large radii. For physical reasons discussed in Sect. 3.1, the focus of the paper is on the first two families. 



For simplicity, we consider two families of distribution func- 
tions. The first family is defined as 

f* T (I) = Ae- aE ° [e-«'- E °> - 1 + a(I - £„)] (21) 

if E < Eq and f^ T (J) = otherwise, so that both f^ T (J) and 
its derivative with respect to E are continuous. We refer to this 
truncation prescription as Wilson truncation (hence the subscript 
WT) because, in the limit of vanishing internal rotation (to — > 0), 
this family reduces to the spherical limit of the distribution func- 
tion proposed by Wilson ( 119751 ). The superscript d in Eq. d2"TT > 
indicates the presence of differential rotation. 

The second family is defined by the distribution function 

// r (/) = Ae -°Eo e -a(i-E ) (22) 

if E < £0 and fp T (J) = otherwise; therefore the function, 
characterized by plain truncation (hence the subscript PT), is 
discontinuous with respect to E. In the limit of vanishing in- 
ternal rotation, it reduces to the spherical limit of the function 
proposed by Prendergast & Tomer (1970), which leads to the 
truncated isothermal sphere (see also Woolley & Dickens [l962l ). 
A summary of the main properties and definitions of the rele- 
vant nonrotating limit of the two families of models is presented 
in AppendixlAl 

In both cases the distribution functions are positive definite 
/wt-C-Oi fpjiD ^ 0, by construction. Curiously, a naive extension 
of King models / = ft(J) with a similar truncation in energy 
alone would define a distribution function that is not positive 
definite in the whole domain of definition. 

Sharp gradients or discontinuities in phase space (such as the 
ones associated with the truncation prescription of fp T (T)) are 
expected to be associated with evolutionary processes dictated 
either by collective modes or by any small amount of collisional- 
ity. Therefore, the first truncation prescription, corresponding to 
a smoother distribution in phase space, is to be preferred from a 
physical point of view as the basis for a realistic equilibrium con- 
figuration (in principle, we might have referred to even smoother 
functions; see Davoust fl977l ). In addition, a full analysis of the 
configurations defined by f^ T (J) shows that this family of mod- 
els exhibits a number of interesting intrinsic and projected prop- 
erties, more appropriate for application to globular clusters, with 
respect to the models defined by fiJT). 

Therefore, the following Sects. [4]and[5]are devoted to the full 
characterization of the family of models defined by f£ T (T) (for 



a summary of the properties of the families of models studied 
in the present paper, see Table 1). The intrinsic properties of the 
family of models defined by fp T (I) are summarized in Appendix 
B. In this investigation, we decided to briefly mention and to 
keep also the second family not only because it extends a well- 
known family of models, but also because it allows us to check 
directly an important aspect of model construction that had been 
noted by Hunter d 19771 ). This is that the truncation prescription 
affects the density distribution in the outer parts of the models 
significantly. This general point is even more important if we are 
interesting in modeling the outermost regions of globular clus- 
ters (and the issues that are often discussed under the name of 
studies of "extra-tidal light", e.g., see Jordi & Grebel [2010l ). 

3.2. The construction of the models 

The construction of the models requires the integration of the rel- 
evant Poisson equation, supplemented by a set of boundary con- 
ditions equivalent to the one described in Sect. 2 for rigidly ro- 
tating models. In this case, we obtain the solution by means of an 
iterative approach, based on the method proposed by Prendergast 
& Tomer d 19701 ), in which an improved solution iff {n) of the 
dimensionless Poisson equation is obtained by evaluating the 
source term on the right-hand side with the solution from the im- 
mediately previous step (for an application of the same method 
to the construction of configurations shaped by an external tidal 
field, see Sect. 5.2 in Paper I) 

V 2 (//"> = — -d(r,e,^"- l) ) ; (23) 
Po v ; 

here the dimensionless escape energy is 

<A(r) = a[E - 0> c (r)] , (24) 

the dimensionless radius is defined as r = r/ro, with the same 
scale length introduced in Eq. ©, and po indicates the dimen- 
sionless central density. The relevant density profile p = Ap, 
with A defined as in Eq. (0, results from the integration in ve- 
locity space of the distributions function defined by f^ T (I). It 
is clear that the general strategy for the construction of the self- 
consistent solution is applicable also to the density derived from 
fp T (J). The iteration is seeded by the corresponding spherical 
models, that is, the Wilson and the Prendergast-Tomer spheres 
respectively, and is stopped when numerical convergence is 
reached (see Appendix ICl for details). At each iteration step, the 
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Fig. 9. Two-dimensional parameter space, given by central ro- 
tation strength x vs - concentration of differentially rotating 
models defined by f£ T (T); the remaining parameters are fixed 
at b = c = 1. The upper solid line marks the maximum ad- 
mitted values of x, for given values of concentration "P, that is 
the underlying area indicates pairs Q¥,x) for which models can 
be constructed. The intermediate and the lower solid lines mark 
the values of <2)/ci> max = 0.4,0.2, respectively. The gray, wide- 
striped, and thin-striped areas represent the extreme, rapid, and 
moderate rotation regimes, respectively. 




Fig. 10. Values of the ratio between ordered kinetic energy and 
gravitational energy t = K ori t/\W\ for selected sequences of dif- 
ferentially rotating models characterized by *T = 1, 2, .., 7 (from 
bottom to top) and x m the range [Q,Xmax]\ the remaning pa- 
rameters are fixed atb-c- I. The arrows mark the threshold 
values of x for the moderate and rapid rotation regimes, illus- 
trated in Fig. [9] For comparison, the dashed line represents the 
values of t for the sequence of Maclaurin oblate spheroids (with 
e < 0.92995; this eccentricity value corresponds to spheroids 
with maximum value of the rotation parameter Xmax = 0.1 1233). 



scheme requires the expansion in Legendre series of the density 
and the potential 



(25) 



(26) 



i/K (r) are therefore 

— 2 d 1(1+1) 

dr 2 r dr f 2 



The associated Cauchy problems for the radial functions 



, (hi 9 „(„_!) 

n = —rp] 

Po 



supplemented by the following boundary conditions 
(^"'(O) = *fV2 , 
1/^(0) = 0, for/ *0 
^'(0) = *f(0) = , 



(27) 

(28) 
(29) 
(30) 



where *P is the depth of the dimensionless potential well at the 
center. By using the method of variation of arbitrary constants, 
the radial functions can be expressed in integral form as follows 

t// ( "\r) = *PV2- — P r'ff?- l \r)df> 
Po [Jo 

4 f r' 2 p "-\r')dr'} , (31) 
r Jo 



Xoo 
r' l - l pf-'\r')dr' 



Vl (21+ l)p 

^ fr'^pf-Xr'W' 
r Jo 



(32) 



The factor V2 appearing in Eqs. (l28l i and (|3TT > is due to the nor- 
malization assumed for the Legendre polynomials. 

3.3. The parameter space 

Much like in the case of rigidly rotating models, in both fami- 
lies fw T (I) and fp T (I), the resulting models are characterized by 
two scales, associated with the positive constants A and a, and 
two dimensionless parameters (*T,^f), measuring concentration 
and central rotation strength, respectively. In addition, two new 
dimensionless parameters, namely c and 



brl c a c 



(33) 



determine the shape of the rotation profile. Variations in the pa- 
rameter b and c are found to be less important. Minor changes in 
the model properties are found up to b, c » 4, above which the 
precise value of c has only little impact on the properties of the 
configurations. 

For each family of differentially rotating models, for given 
values (*T, b, c) there exists a maximum value of the central ro- 
tation strength parameter Xmax, corresponding to the last config- 
uration for which the iteration described in Sect. l3.2l converges. 
Such maximally rotating configurations exhibit highly deformed 
morphologies, characterized by the presence of a sizable central 
toroidal structure, which will be described in detail in the fol- 
lowing sections. 
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As for the parameter space of rigidly rotating models, it is 
useful to introduce different rotation regimes, defined on the ba- 
sis of the deviations from spherical symmetry introduced by the 
presence of differential rotation. With particular reference to the 
parameter space of the models defined by f^ T (I), we introduce 
some threshold values in central dimensionless angular velocity, 
which, in this family of models, is related to the central rota- 
tion strength parameter^ by the relation a> = 3^'^ 2 , as for the 
rigidly rotating models. In particular, configurations in the mod- 
erate rotation regime have u>/cj max < 0.2 (the thin-striped area in 
Fig.|9]l, are quasi-spherical in the outer parts, while they are pro- 
gressively more flattened when approaching the central region, 
as the value of x increases. For the models falling in this rotation 
regime, the central toroidal structure is absent or, when low val- 
ues of the concentration parameter *P are considered, not signifi- 
cant. Configurations with 0.2 < u>/cj max < 0.4 (the wide-striped 
area in Fig. |9]l are defined as rapidly rotating models. The ex- 
treme rotation regime is defined by the condition a>/a> max > 0.4 
(the gray area in Fig. [5); in this case, the models always show 
a central toroidal structure, which becomes more extended as 
the central rotation strength increases. In particular, in the last 
regime, the entire volume of a configuration is dominated by the 
central toroidal structure. 

As for the rigidly rotating models described in Sect. 2, a 
global kinematical characterization is offered by the parameter 
t = K or dl\W\, defined as the ratio between ordered kinetic en- 
ergy and gravitational energy. Figure [10] illustrates the relation 
between the two parameters, for models with selected values 
of X Y, b, and c. Note that the transition from rapid to extreme 
rotation corresponds to values of the parameter t in the range 
[0.075,0.135] (the precise value depends on the value of x ¥), 
A naive application of the Ostriker & Peebles ( 119731 ) criterion, 
which states that axisymmetric stellar systems with t > 0.14 are 
dynamically unstable with respect to bar modes, would suggest 
that the majority of the models in the extreme rotation regime are 
dynamically unstable. A detailed stability analysis of the con- 
figurations in the three rotation regimes has been performed by 
means of specifically designed N-body simulations and will be 
presented in a separated paper (Varri et al. 2011, Varri et al. in 
preparation). 

To provide a systematic description of the intrinsic and pro- 
jected properties, we will study the equilibrium configurations 
as sequences of models characterized by a given value of con- 
centration T*, in the range [1,7], and increasing values of^-, up 
to the maximum value allowed; such sequences are constructed 
by fixing b = c = 1 (unless otherwise stated). 

4. Intrinsic properties of the differentially rotating 
models 

4.1. The intrinsic density profile 



where the function in the integrand is defined as 



g(s, t, r, 0) = exp 



5 to r sinty 
1 +b[ V2sfrsin6>] 2e J 



(35) 



for completeness, we note that the two dimensionless variables 
in the double integral can be expressed in terms of the previous 
variables as t — cos ft and s = av 2 /2. Because the distribution 
function depends only on the energy E and the z-component of 
the angular momentum J z , the resulting models are axisymmet- 
ric and therefore the density profile depends only on the radius f 
and the polar angle 9. The density depends on the spatial coordi- 
nates explicitly and implicitly, through the dimensionless escape 
energy if/(f, 8); such explicit dependence is the reason why, in 
this case, the isodensity, equipotential, and isobaric surfaces are 
not in one-to-one correspondence, at variance with the family of 
rigidly rotating models. The presence in Eq. (l34l of the terms 
with fractional powers of i[> is due to the adopted truncation in 
phase space; in particular, it is directly related to the presence of 
the terms e~ aE ° [- 1 + a(I - E )] in Eq. (TUT) . 

The central value of the density profile depends only 
on the concentration parameter T* and is given by pwrfl = 
(2/5)e y(7/2, V P), consistent with the central value of the den- 
sity profile obtained in the nonrotating limit piyr,s( v P) (see 
Eq. ( IA.4I ) in Appendix [A}. This result corresponds to the fact 
that the integral I(E, J z ) (see Eq. ( l20b ) reduces to the Jacobi in- 
tegral for small values of J z , which implies that the rotation is 
approximately rigid in the central regions of a configuration (see 
the next subsection for details); therefore, the mean rotation ve- 
locity vanishes at the origin and the density distribution reduces 
to its nonrotating limit. 

The double integral in Eq. ( l34b requires a numerical inte- 
gration (see Appendix ICl for details). Some insight into the be- 
havior of the density profile can be gained by calculating the 
relevant asymptotic expansions in the central region and in the 
outer parts. Around the origin, up to second order in radius, the 
density profile reduces to 

p WT (r, e, *P) = pwr,o + \^y{^> T ) [ 9 * si " 2 6 

d 2 iff 



df 2 



r 2 +0(r 4 ), 



(36) 



which depends explicitly on the concentration *T and the cen- 
tral rotation strength^, and implicitly on the parameters b and c, 
through the second order derivative of the escape energy evalu- 
ated at f = 0. Such derivative can be calculated from the radial 
functions given in Eqs. (|3TI)-(|321) 



df 2 



The relevant density profile is obtained by integration in velocity C2/5\^ 2 / 2 \ 

space of the distribution function f^ T {I) (see Eq. (UTTi). It is con- =_3 + ~2~\2j ( _1+ 3cos0j, (37) 



venient to introduce in the velocity space a spherical coordinate 
system (v, fi, A), in which v is magnitude of the velocity vector, 
while fi and A are the polar and azimuthal angle, respectively. 
After some manipulation, the density profile can be expressed in 
dimensionless form as 



where the quantity 

is r + °° i 

5po Jo r> 



(38) 



pwr{f,G,ij/) = 



-e r 
4 



dse 



•s.l/2 



J dt g(s, t, r, 6) 



^3/2 _ 



(34) 



depends implicitly on the parameter ^ through the function p2(r) 
and po. The quantity C 2 is negative-definite since the quadrupole 
radial function of the density piif) is negative on the entire do- 
main of definition of the solution of the Poisson equation; the 
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Fig. 11. Intrinsic density profiles (normalized to the central 
value) evaluated on the equatorial plane (solid lines) and along 
the S-axis (dashed lines) for a sequence of differentially ro- 
tating models defined by f^, T (I), with *F = 2 and x - 
0.04,0.36,1.00,1.96,3.24 (from right to left; slower rotating 
models are more extended); the remaining parameters are fixed 
atb - c - I. 



sign of the quadrupolar function is negative because the config- 
urations in our family of models are always oblate. In passing, 
we also note that the the expansion around r — of the escape 
energy up to second order in radius is given by 



ip(f, 0) = *P + 



1 <9V 

2 ~dP 



r 2 + 0(f) . 



(39) 



Since the boundary of a configuration is defined by the con- 
dition if/(f, 0) = 0, the density profile in the outer parts can be 
evaluated by performing an expansion with respect to if/ <s 1 



pwrif, 8, iff) 



zxr 1 sin 2 1 



if^ 2 +0(if/" z ); 



7/2x 



the terms with fractional powers of if/ that appear in Eq. d34i > 
cancel out with the first terms of the expansion of the double 
integral (which reduces to the incomplete gamma function). 

The density profiles evaluated on the principal axes for a se- 
quence of models with *P = 2, b - c - 1, and increasing values 
of the rotation strength parameter x are illustrated in Fig. QT] 
the corresponding dimensionless escape energy profiles are dis- 
played in Fig. [12] Configurations characterized by moderate ro- 
tation show monotonically decreasing profiles, whereas models 
with rapid rotation have the maximum value of the density pro- 
file in a position displaced with respect to the origin; in the ex- 
treme rotation regime, also the maximum value of the escape en- 
ergy is off-centered. The sections of the isodensity and equipo- 
tential surfaces (presented in the first two rows of Fig. [TBI clearly 
show that the offset of the density peak corresponds to the exis- 
tence of a curious toroidal structure; the condition for the exis- 
tence of such structure is discussed in Sect. 14.31 

4.2. The intrinsic kinematics 

The calculation of the first order moment in velocity space of 
the distribution function f£ T (T) confirms that only the azimuthal 



Fig. 12. Dimensionless escape energy (normalized to the central 
value) evaluated on the equatorial plane (solid lines) and along 
the z-axis (dashed lines) for the sequence of differentially rotat- 
ing models displayed in Fig.fTTI 



component of the mean velocity is nonvanishing. The mean ve- 
locity in dimensionless form 



{v<t,) WT (r,6,iff) = 



-In g(s,t, r,6)] 



dss J dtt [g(s,t,r,6)e~ s e 



(41) 



can be calculated by numerically. As for the density profile, it 
is useful to evaluate the asymptotic expansion of Eq. (BTI) in the 
central regions and in the outer part of a given configuration. By 
performing a first order expansion in the radius with respect to 
the origin, we found that the mean rotation velocity reduces to 
the following expression 



(40) (%) WT (f, 0) = 3x 1/2 r sin + O^) 



(42) 



which corresponds to rigid rotation, with dimensionless angular 
velocity oj = 3x 1 ^ 2 ', the expression does not depend on the con- 
centration parameter, consistent with the asymptotic properties 
of the integral I(E, J z ). 

In the outer parts of the models the mean rotation velocity 
is of order 0{if/) and thus vanishes at the boundary. The relevant 
numerical coefficients depend on the value of the parameter c; 
for example, for c = 1 we have 



2(2 + 6br 2 sm 2 + 9xr 2 sm 2 0) 
(v^) WT (r,e,if/) = — , if/ 



(43) 



The mean rotation velocity profiles on the equatorial plane 
for a selected sequence of models are displayed in Fig. |T4j as 
the value of the central rotation strength parameter increases, the 
mean velocity profile becomes steeper in the inner parts and the 
maximum value increases; the central part of the profiles is well 
approximated by rigid rotation, as in Eq. d42l) . 

By evaluating the second-order moments in velocity space 
of the distribution function f£ T (T), it can be easily shown that, 
in our coordinate system, the pressure tensor p,j - (A/a)pij is 
diagonal and that the radial component is equal to the polar one 
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Fig. 13. Meridional sections of the isodensity, equipotential, isovelocity L and isobaric surfaces (from top row to bottom row) for a 
sequence of four differentially rotating models characterized by *F = 2, b — c — 1, and^f = 0.04, 0.16, 0.36, 1.0 (from left column 
to right column; note the change in the scale of the axes). The first and second models are in the moderate rotation regime, the third 
has rapid rotation, and the last represents the beginning of the extreme rotation regime. 
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Fig. 14. Mean rotation velocity profiles on the equatorial plane 
(solid lines) for the sequence of differentially rotating models de- 
fined by f* T (l), with ¥ = 2 and^ = 0.04, 0.36, 1.00, 1.96, 3.24 
(the same sequence displayed in Figs. fTTI and fT2l from bottom 
to top). Dashed lines indicated the asymptotic behavior in the 
central regions, which corresponds to rigid rotation (see Eq. (l42l 
and Sect. 1431), 



Pn = Pee- Therefore, in the following, only the nontrivial com- 
ponents are discussed. The radial and azimuthal components in 
dimensionless form are given by 

3 r +1 

p Wirr (f,6,f) = -z* dseT s s 3/2 dt{\-t 2 )g(s,t,r,0) 

4 Jo J-i 



.^5/2 _ 1^7/2 

5^ 35^ 



(44) 



Pw44,(r, 6, if/) 



dse 



-s c 3/2 



J dt t 2 g(s, t, r, 6) - 



J<A 5/2 -^<A 7/2 -Pwr<V0>wr 



(45) 



where the presence of the terms with fractional powers of if/ 
should be interpreted as in Eq. (l34b . The expansion up to sec- 
ond order in radius gives 

p %rr {f, 8, ¥) = pwt,o + p*y \1, [9 X sin 2 
r 2 + 0(r 4 ) . 



(46) 



We also found that, to second order in radius, the pressure 
is isotropic pw,^(f, 0, *P) = Pw,rr(?, W). The components of 
the pressure tensor evaluated at the origin reduces to pwr,o - 
(4/35)e' i V(9/2, v f'), consistent with the value obtained in the 
nonrotating limit Pwt,s(^) (see Eq. ( IA.5t in AppendixlAl. 

The asymptotic behavior of the pressure tensor components 
in the outer part can be written as 



18 

p %rr (t 6, if/) = — X r 2 sin 2 6if/ 7/2 + 0(if/ 9/2 ) , 
PwM" r > > <A) = 5 ^Xr 2 sin 2 6i(/ V2 + 0(tf/ 9/2 ) , 



(47) 
(48) 




<b 0.2 



Fig. 15. Squared velocity dispersion profiles for the azimuthal 
(solid lines) and radial component (dashed lines) of the sequence 
of differentially rotating models displayed in Fig.[[4](from right 
to left). The profiles are evaluated on the equatorial plane. 



-0.5 



T 




-tj 



Fig. 16. Radial profiles of the anisotropy parameter for the se- 
quence of differentially rotating models displayed in Fig. Q3] 
evaluated on the equatorial plane (from right to left). All models 
are characterized by a — (isotropy) at the center and a = —2 
(tangentially-biased anisotropy) at the boundary. 



respectively. 

The relation between the dimensionless pressure and veloc- 
ity dispersion tensor is given by <x 2 = pij/p. The profiles of the 
radial and azimuthal component of the velocity dispersion ten- 
sor of a selected sequence of models, evaluated on the equatorial 
plane, are displayed in Fig. Q3] For configurations in the mod- 
erate rotation regime the profiles are monotonically decreasing 
with the radius, with some variations in the slope in the interme- 
diate and external parts, while for configurations in the rapid and 
extreme rotation regimes, their peak is off-centered. The sections 
in the meridional plane of the isobaric surfaces, defined with re- 
spect to the trace of the pressure tensor, show the presence of a 
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Fig. 17. Radial profiles of the anisotropy parameter for selected sequences of differentially rotating models, evaluated on the equa- 
torial plane. Top panels: the sequences are characterized by *P = 2, x — 0.04,0.16,0.36,0.64, 1.00 (in each panel, from right to 
left); the left and right panels show the effect of increasing of parameters c and b, respectively. Bottom panels: the sequences are 
characterized by b — c — 1, x — 0.04,0.16,0.36,0.64, 1.00 (in each panel, from right to left); the left and right panels show the 
effect of decreasing and increasing the concentration parameter *P, respectively. 



central toroidal structure for the fast rotating models of the se- 
quence (see the last row in Fig. [TBI. 

The intrinsic kinematics can be further characterized by 
means of the anisotropy parameteiQ , defined as 



1 



P<M> 

Prr 



1 



o~2 



(49) 



From the asymptotic behavior of the pressure tensor components 
in the central region (see Eq. d46*il). we find that a — > for r — > 0, 
while from the expansion in outer parts (see Eqs. d47l)-(l4"8ii). we 
find that a — > -2 as we approach the boundary. These limit- 
ing values do not depend on the dimensionless parameters that 
characterize the family f^ T (I). In other words, the central re- 
gion of the configurations is always characterized by isotropy 
in velocity space, while the regions next to the boundary show 
a strong tangentially-biased pressure anisotropy. The radial pro- 
files of the anisotropy parameter for a selected sequence of mod- 
els are displayed in Fig. [16] The values of a in the interme- 
diate region of a given configuration depend on the values of 
the relevant dimensionless parameters Q¥, b, c). In fact, from an 
exploration of the entire four-dimensional parameter space, we 
found that, by increasing the value of c, the portion of a model 
dominated by radial anisotropy becomes slightly less extended, 
while, if the value of b is increased, it increases (see Fig. [P71 



4 In the literature, the anisotropy parameter is often defined as /? = 
1 - (trL + a 2 eg )l2a 2 r ; for axisymmeric systems, for which a 2 m = cr 2 r , 
the relation with the parameter adopted in the present paper is given by 
a = 2/3. Here we prefer to keep the same notation used in van Albada 
119821 Zocchi et al. 2012 and in other articles. 



top panels). In turn, by increasing the value of the concentration, 
as measured by y ¥, the models are characterized by a significant 
radially-biased pressure anisotropy, which appears also in the in- 
termediate region. By decreasing the value of the concentration, 
the intermediate region turns out to be dominated by tangential 
anisotropy, like in the outer parts (Fig. [17] bottom panels). 



4.3. The condition for the existence of the central toroidal 
structure 

In general, configurations with rapid rotation may exhibit highly 
deformed morphologies. In the context of rotating fluids, the 
regime of strong differential rotation has been successfully ex- 
plored, at least for polytropes. Stoeckly d 19651 1 and Geroyannis 
( 1 19901 ) found that rapidly rotating polytropic fluids show a cen- 
tral toroidal structure; in addition, a self-consistent method for 
the construction of rapidly differentially rotating fluid systems 
with a great variety of shapes has been proposed by Hachisu 
( 1986). In stellar dynamics, this regime has been rarely explored: 
Lynden-Bell (119621 ) and Prendergast & Tomer d 19701 ) noted that 
some models with strong differential rotation show the density 
peak in a ring on their plane of symmetry. 

The existence of a central toroidal structure in a given model 
of our family can be studied by using the asymptotic expansion 
of the density in the central regions, expressed in Eq. ( f36b . In 
fact, if the second order derivative of the density with respect to 
the radius is positive, then the maximum value of the radial den- 
sity profile is displaced from the geometric center of the config- 
uration, so that a central toroidal structure is formed. Therefore, 
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Fig. 18. Two-dimensional parameter space of differentially rotat- 
ing models C¥,x)', the remaining parameters are fixed at values 
b = c = 1 . The upper solid line marks the maximum admitted 
values of x, as in Fig. [9] The lower solid line marks the values 
of x at which the models start to exhibit a central toroidal struc- 
ture; the dashed line marks the upper limit of the condition for 
the existence of the central toroidal structure (see Eq. (BTTl). 



the relevant condition on the sign of the derivative can be ex- 
pressed as 



9*-sin 2 6»-3 + y (-1 + 3 cos 2 0) > , 



(50) 



which, on the equatorial plane (i.e., at 6 — n/2), reduces to a 
simple condition for the central rotation strength parameter 



1 Co (5 



y > — I 

A 3 18 \ 2 



1/2 



(51) 



The expression on the right-hand side depends implicitly on the 
dimensionless parameters T*, b, and c through the quantity C2, 
defined by Eq. d38T >. Since C2 is negative-definite, the condition 
X > 1/3 provides the upper limit to Eq. (15 It . Actually, the mod- 
els characterized by values of x immediately above the thresh- 
old given by Eq. ( BTT l, show a very small central toroidal struc- 
ture, with a shallow increase of the density with respect to the 
geometric center of the configuration (e.g., for the model with 
¥ = 2 and x = 0.36, illustrated in Fig. QT| the density peak, at 
the center of the toroidal structure, is merely log(p/po) = 0.014 
and the central structure itself is barely visible in the meridional 
sections of the isodensity surfaces, depicted in the third panel of 
the first row of Fig.fTJll. To construct a model with a sizable cen- 
tral toroidal structure, the value of the parameter x should be at 
least twice the threshold given by the above-mentioned condition 
(e.g., see the last panel of the first row of Fig.Qj). In conclusion, 
by comparing Figs. [9] and [18] we note that, for configurations 
in the moderate rotation regime (i.e., with u>/u> max S= 0.2), the 
central toroidal structure is absent or very small, while starting 
from the rapid rotation regime, the structure is always present 
and becomes more extended as the value of the central rotation 
strength increases. 

The presence of the central toroidal structure can be also be 
interpreted in terms of to the intrinsic kinematical properties the 




0.30 



Fig. 19. Squared angular velocity (solid lines) evaluated in the 
central parts of the equatorial plane for the first three differen- 
tially rotating models displayed in Fig. [13] that is with = 2, 
b = c = 1, and^- = 0.04, 0.16, 0.36 (from bottom to top). Dashed 
lines indicate the asymptotic behavior at small radii, character- 
ized by solid-body rotation, that is constant angular velocity. 
Dotted lines represent the angular velocity associated with the 
circular orbit of a single star, evaluated on the equatorial plane 
for the same models (from top to bottom). For the first two mod- 
els Gj c > at, while for the third cb > cb c in the inner part, where 
the toroidal structure is present (the black filled circle marks the 
position where a> = u> c ) . 

models. From the radial component of the Jeans equation, ex- 
pressed in dimensionless spherical coordinates, 



p dr r 



1 dip <v>) 2 



dr 



IdSrl 



0-7., 



dr 



(52) 



the sign of the first order derivative of the density with respect to 
the radius (on the left-hand side of Eq. (l52l ) depends on (i) the 
difference between the angular velocity associated with the cir- 
cular orbit of a single star Gj c = [-(l/r)d- r ip] 1 ^ 2 and the angular 
velocity of the model <2>, associated to the mean rotation veloc- 
ity (v^) = <br sin 6, (ii) a more complex pressure term (in square 
brackets on the right-hand side of Eq. (f52l). To check for the 
presence of the central toroidal structure, it is sufficient to study 
the Jeans equation in the central region of the model. Therefore, 
by inserting the relevant asymptotic expansions for the mean ve- 
locity and the escape energy (see Eqs. d42l) and d39i>). the term 
(i) reduces to the expression indicated on the left-hand side of 
Eq. ( l50b . and, by inserting the relevant asymptotic expansions 
for the pressure tensor components, the term (ii) can be written 
as 

1/2 



9;fsin 2 6»-3 + y (~j (-l + 3cos 2 6») 



5 7(9/2,^)7(5 /2,¥) 
7 r(7/2, >F) 2 



(53) 



By combining the asymptotic expressions of terms (i) and (ii), 
it follow^] that the requirement of a positive density gradient 

5 The coefficient given by 1 - (5/7)[y(9/2, v P)y(5/2, v F)]/y(7/2, Y) 2 
is nonnegative for every value of *P. 
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Fig. 20. Contour maps of the surface density, mean line-of-sight velocity, and line-of-sight velocity dispersion (from top row to 
bottom row) of the differentially rotating models with *F = 2, b = c = 1, and^- = 0.04,0.16,0.36, 1.00 (from left to right, as in 
Fig- US, projected along the i-axis of the intrinsic coordinate system ("edge-on" view, so that jtp and yp correspond to the y and z 
axes of the intrinsic system, respectively). In the panels in the first row, solid lines represent the isophotes corresponding to selected 
values of S/Sq in the range [1.02, 10~ 7 ]; only the last (fastest rotating) model shows values E/Eo > 1, when the toroidal structure 
appears. Panels in the second row illustrate the contours of the dimensionless rotation velocity in the range [0.5, 10~ 5 ] at intervals 
of 0.05 (from left to right, the values of the innermost contours are 0.25, 0.35, 0.4, 0.45, respectively.). Panels in the last row show 
the contours of the projected velocity dispersion in the range [0.4, 10~ 5 ] at intervals of 0.05 (the values of the innermost contours 
are 0.3 for the first three models and 0.4 for the last one). 



on the left-hand side of the Jeans equation is equivalent to the 
condition expressed by Eq. (1501 . 

In conclusion, we have independently tested the validity of 
the condition for the existence of the central toroidal structure 
derived at the beginning of this section. We have also shown that 
the requirement of positivity of term (i) in the radial Jeans equa- 
tion is a necessary and sufficient condition for the presence of 
the central toroidal structure. In other words, the central toroidal 
structure exists if and only if, in the central region, the angu- 
lar velocity associated with the internal rotation is higher than 
the angular velocity associated with the circular orbit of a single 
star. Figure [19] shows the relevant angular velocities for the first 
three models of the sequence considered in Fig.Q~3j it is apparent 
that, for the configuration in which the central toroidal structure 
is present, u > a> c . This result strictly depends on the adopted 
truncation in phase space for the distribution function f£ T (T); in 
fact, in Appendix IE1 we show that for the alternative distribution 
function fp T (I) , the condition on the angular velocities is a nec- 



essary but not sufficient condition for the existence of the central 
toroidal structure. 

4.4. The condition for maximally rotating models 

If we consider a sequence of models with fixed values of V F, b, c 
and increasing values ofx, the central toroidal structure becomes 
progressively more extended and characterized by a larger aper- 
ture angle (i.e., the angle spanned at small radii by the isodensity 
surface that represents the boundary of the central toroidal struc- 
ture). For high values of the central rotation strength parameter 
X, a smaller central toroidal structure appears also in the equipo- 
tential surfaces (see the escape energy profile of the fastest ro- 
tating model in Fig. [12}. These two structures can be charac- 
terized in terms of P and 8 C {, defined as the complements of the 
semi-aperture angle of the inner cusp of toroidal structures in the 
equipotential and isodensity surface, respectively. In fact, since 
the boundary of the central cusp of the toroidal structures, is de- 
fined by if/(f, 6) = ^P, and p(r, 6, i/S) - po, by using the relevant 
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Fig. 21. Surface density profiles (normalized to the central value) 
of the differentially rotating models with ¥ = 2, b - c - 1, 
and^ = 0.04, 0.16, 0.36, 1.00 (from right to left; same sequence 
illustrated in Fig. l20b . evaluated along the xp-axis (solid lines) 
and yp-axis (dashed lines) of the projection plane ("edge-on" 
view). 



Fig. 22. Ellipticity profiles, as functions of the semi-major axis 
of the projected image ap, of the first three differentially rotat- 
ing models illustrated in Fig. [2T](from bottom to top; "edge-on" 
view). Dots correspond to the isophotes shown in the first row of 
Fig.[20]and arrows mark the position of the half-light isophote. 



asymptotic expansion up to second order in f given in Eqs. ([39l > 
and d36l ). the following expressions for the angles are obtained: 



cos 2 9 n 



3 + (5/2)'/ 2 C 2 /2 
(3/2)(5/2) 1 / 2 C 2 



9 X 



2 3 + (5/2)'/ 2 C 2 /2. 

cos Gd = ttz 

(3/2)(5/2)^C 2 -9 X 

The two angles are related, because 



(3/2)(5/2) l / 2 C 2 cos 2 9 p -9 X 

cos 9 d - 



(54) 



(55) 



(3/2)(5/2)i/2c 2 - 9* 



(56) 



Note that 9 p and 9 d decrease as the values of \C 2 \ and^- increase 
(we recall that C 2 is negative-definite), that is the corresponding 
semi-aperture angles become larger. 

For high values of the central rotation strength X , which im- 
ply a high degree of quadrupolar deformation, as measured by 
the quantity |C 2 |, it is easy to see that cos 2 6^ — > 1/3; interest- 
ingly, we also found (numerically) that a limiting value exists 
also for 9 d , given by cos 2 9 d — > 2/3. By inserting the limiting 
value of 9 P in Eq. d56t . the condition for the maximum value of 
the angle 9d (i.e., cos 2 9 d < 2/3) can be translated in a simple 
condition for the central rotation strength parameter 



X< 



6 2; 



1/2 



\c 2 



(57) 



which basically provides a limit to the deformation induced 
by the central rotation itself. The previous condition has been 
checked numerically and we found that, for values of x above 
the threshold, the iteration scheme used for the solution of the 
Poisson equation does not converge. 



5. Projected properties of the differentially rotating 
models 

5.1. The surface density profile 

The calculation of the projected properties has been performed 
by following the same projection rules described in Sect. I2.4I 
that is the line of sight corresponds to the £p-axis of a new frame 
of reference in which the projection plane is denoted by (xp, yp). 
In particular, we studied in detail the "edge-on" view, that is the 
projection along the x-axis of the intrinsic frame of reference 
(zp = x, Xp = y, and yp = z, i.e., the viewing angles of the rota- 
tion matrix are 9 — n/2, <p = 0). The dimensionless surface den- 
sity distribution £(ip, yp) has been calculated by numerical in- 
tegration of Eq. (fT9] l on an equally-spaced square cartesian grid 
on the projection plane. 

The first row of Fig. [20] shows the contour maps of the sur- 
face density of selected differentially rotating models in the mod- 
erate and rapid rotation regime. As result of projection, the dim- 
ples on the rotation axis, which are prominent in the correspond- 
ing intrinsic isodensity surfaces (see Fig. [13] for the meridional 
sections of the same sequence of models), are less pronounced. 
In addition, the central toroidal structure in the projected density 
distribution is visible only if it has a reasonable size in the intrin- 
sic density distribution (see third and last model of the sequence 
illustrated in Fig.l20l. 

The surface density profiles of the same sequence of differ- 
entially rotating models, evaluated along the principal axes of 
the projection plane, are presented in Fig. [2T| For the config- 
urations in which the central toroidal structure in projection is 
absent, we calculated the relevant ellipticity profiles, as func- 
tions of the semi-major axis of the projected image ap, and they 
are reported in Fig. [22] As expected, the configurations in the 
moderate rotation regime are characterized by nonmonotonic el- 
lipticity profiles, while models in the rapid rotation regime have 
monotonically decreasing profiles; to some extent, this morpho- 
logical feature is complementary to that of the uniformly rotating 
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Fig. 23. Mean line-of-sight rotation velocity profiles of the dif- 
ferentially rotating models with *F = 2, b = c = 1, and 
X = 0.04,0.16,0.36, 1.00 (slower rotating models are more ex- 
tended; same sequence illustrated in Fig. I2TI 1. evaluated along 
the xp-axis of the projection plane ("edge-on" view). 

models, in which the configurations are always characterized by 
monotonically increasing ellipticity profiles. Interestingly, the 
behavior of the ellipticity profiles does not necessarily correlate 
with the mean line-of-sight velocity profiles (see next subsec- 
tion), that is configurations with a nonmonotonic velocity profile 
may have a monotonic ellipticity profile. 

In addition, the isophotes of models show clear departures 
from a pure ellipse, which, at variance with the family of rigidly 
rotating models, can be characterized as a "boxy" overall trend, 
particularly evident in the intermediate parts of the configura- 
tions (see third and fourth panels of the first row of Fig. 1201. 

5.2. The line-of-sight velocity distribution 

The projected velocity moments are calculated by integrating 
along the line of sight (weighted by the intrinsic density) the 
corresponding intrinsic quantities. As for the surface density dis- 
tribution, we studied in detail the kinematics projected along the 
x-axis of the intrinsic coordinate system. Therefore, the dimen- 
sionless mean line-of-sight velocity and the line-of-sight veloc- 
ity dispersion can be written as 



(v P }(xp,y P ) 
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<v>> x P 
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Z P + X 



Xp 



L) 1/2 



j(z 2 P & 2 rr + 



(58) 



(e+4) 1/2 



(59) 



Selected mean line-of-sight velocity profiles (for the se- 
quence of models presented in Fig. 12Tb. evaluated along the xp- 
axis of the projection plane, are illustrated in Fig. [23] where the 
sign of the mean velocity in the two half-planes is consistent with 
Eq. (l58l . As the value of the central rotation strength parame- 
ter increases, the slope of the inner part of the profile becomes 



CTp 




rp 

Fig. 24. Squared line-of-sight velocity dispersion profiles of the 
differentially rotating models illustrated in Fig. [23] evaluated 
along the ip-axis (solid lines) and yp-axis (dashed lines) of the 
projection plane ("edge-on" view). 




Xp 

Fig. 25. Ratio of mean line-of-sight (rotation) velocity to the 
line-of-sight velocity dispersion of the differentially rotating 
models illustrated in Figs.[23]and[24l evaluated along the xp-axis 
of the projection plane ("edge-one" view). The first and second 
models are in the moderate rotation regime, the third has rapid 
rotation, and the last represents the beginning of the extreme ro- 
tation regime. 



steeper, since the asymptotic behavior of the intrinsic velocity in 
the central regions is approximately that of a rigid rotation, with 
the angular velocity proportional to x 1 ^ 2 (see Eq. d42l); in ad- 
dition, because the entire configuration becomes more compact, 
the radial position of the peak of the velocity profile shrinks pro- 
gressively. 

The line-of-sight velocity dispersion profiles of the same se- 
quence of models, evaluated along the principal axes of the pro- 
jection plane, are presented in Fig. [24] For configurations in the 
moderate rotation regime, the variations in the slope at interme- 
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h h 

Fig. 26. Top panels: histogram (scaled to unity) of the phase space density N(E, J z ) for a differentially rotating model characterized 
by Wilson truncation (left) with *P = 2,x = 0.16, b = c = 1, and for one characterized by plain truncation (right) with *P = 2,x = 
0.36, b — c — 1. Both models are characterized by w/aj mi7A « 0.2. Bottom panels: Lindblad diagrams for the same models presented 
in the top panels. The graphs have been obtained by a Monte Carlo sampling of the relevant distribution functions with 65 536 
particles. 



diate and outer radii, which characterize the azimuthal compo- 
nent of the intrinsic velocity dispersion tensor, are still visible in 
projection, while the inner part of the profile has a flat core. For 
models in the rapid and extreme rotation regime, the maximum 
value of the profile is displaced from the geometric center; in this 
case, the inner part of the profile has a nontrivial gradient. We 
also calculated the ratio of the mean line-of-sight (rotation) ve- 
locity to the line-of-sight velocity dispersion, as a measure of the 
amount of ordered motions compared to random motions (see 
Fig.[25]for relevant profiles evaluated along the ip-axis, for the 
sequence of models discussed above). The models in the moder- 
ate rotation regime show values that are consistent with those ob- 
served in Galactic globular clusters (see Introduction). The con- 
figurations in the rapid and extreme rotation regime are charac- 
terized by higher values of the ratio, even greater than one; such 
high values are measured only in the class of elliptical galaxies 
known as fast rotators (see Davies et al. [1983). Of course, the 
values of this ratio strongly depend on the line of sight on which 
the projection is performed (we recall that here we illustrate the 
results obtained only from the "edge-on" view, which is the most 
favorable for the detection of ordered motions). 



6. Effect of truncation in phase space 

To explore the nontrivial effects of the two options for the trun- 
cation prescription of the family of differentially rotating mod- 
els (see Sect. 13.11 see also Hunter [T977I ), we selected two rep- 
resentative models and studied the relevant phase space density 
N(E, J z ), defined in such a way that the total dimensionless mass 
is given by M — J N(E, J z )dE dJ z , with the specific energy 

E — aE and z-component of the specific angular momentum 
J z - (r x \) z = v^rsind. The histograms of the phase space 
density, constructed by means of a Monte Carlo sampling of the 
distribution functions, are presented in the top panels of Fig. [26] 
We also constructed the corresponding Lindblad diagrams 
(Lindblad 1933), that is the representation in the plane (E,J Z ) 
of orbits in phase space in a given model (see bottom panels of 
Fig. [26l>. For any given J z , the orbit with the lowest energy cor- 
responds to a circular orbit in the equatorial plane; therefore, in 
both families of models, the circular orbits form the lower cuspy 
boundary of the model in the diagram. The upper boundary cor- 
responds to the energy truncation, introduced by the cut-off con- 
stant £■(). Of course, since both families of models are character- 
ized by internal rotation, the distribution of the orbits is asym- 
metric with respect to the z-component of the angular momen- 
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turn, but different truncation prescriptions determine a different 
distribution of the orbits in phase space. In particular, the Wilson 
truncation introduces a sharp depopulation of retrograde orbits 
with intermediate energy, while the plain truncation is associated 
with a smooth decrease of the number of retrograde orbits with 
high energy. Alternative options for the truncation prescription 
in rotating models correspond to a distribution of orbits in differ- 
ent regions of the diagram (see Fig. 1 in Rowley 1988 ). Note that 
a truncation with respect to the Jacobi integral H — E- coJ- only, 
as in our family of rigidly rotating models f' K {H) (see Eq. (Q])), 
corresponds to a region in the diagram bounded by a straight line 
(with co as slope) and the curves corresponding to the energy of 
circular orbits. 

In addition, as the concentration of the models increases, the 
cusp of the lower boundary, given by the energy of circular or- 
bits, becomes sharper; as noted also by Rowley ( 1988 ), this is the 
reason why centrally concentrated configurations with rigid ro- 
tation cannot be very flat, except for the outer parts, as described 
in Sect. 12131 

7. Discussion and conclusions 

In this paper we have constructed two new families of self- 
consistent axisymmetric models of quasi-relaxed stellar systems, 
characterized by the presence of internal rotation (see Table 1); 
a full description in terms of the intrinsic and projected proper- 
ties has been provided. The main results can be summarized as 
follows: 

- Driven by general statistical mechanics considerations, we 
started by constructing a family of rigidly rotating dynam- 
ical models; this family is defined as an extension of the 
King (1966 ) models to the case of axisymmetric equilibria 
flattened by solid-body rotation, with the relevant distribu- 
tion function dependent only on the Jacobi integral. The con- 
figurations have been constructed self-consistently by solv- 
ing the Poisson-Laplace equation for the mean-field poten- 
tial by means of a perturbation method (described in Paper 
I), using a measure of the rotation strength as the expan- 
sion parameter. The two-dimensional parameter space which 
characterizes the family (concentration and rotation strength) 
can be described in terms of two regimes. Models in the 
low-deformation regime are almost indistinguishable from 
the corresponding spherical King models. Highly-deformed 
models are quasi-spherical in the central regions and show 
significant deviations from spherical symmetry in the outer 
parts; in particular, they are flattened toward the equatorial 
plane and exhibit a sort of "disky" appearance. The resulting 
eccentricity profile is a monotonically increasing function of 
radius; the (finite) central value can be expressed analyti- 
cally in terms of the rotation strength parameter. From the 
kinematical point of view, the models are characterized by 
pressure isotropy and cylindrical rotation. 

- In view of possible applications to globular clusters, we have 
constructed a second family of dynamical models, charac- 
terized by differential rotation, designed to be approximately 
rigid in the central regions and to vanish in the outer parts, 
where the imposed energy truncation is effective. In this 
case, the relevant Poisson equation is solved by means of 
a spectral iteration method based on the Legendre expan- 
sion of the density and the potential. The full parameter 
space is now four-dimensional, with two additional param- 
eters, defining the shape of the rotation profile. Three rota- 
tion regimes can be introduced, namely of moderate, rapid, 



and extreme rotation. However, significant variations in the 
structure of the models are primarily associated with con- 
centration and central rotation strength, as for the previous 
family. We explored the properties of the configurations re- 
sulting from two options for the truncation prescription, with 
emphasis on the family which, in the limit of vanishing inter- 
nal rotation, reduces to the spherical limit of the models pro- 
posed by Wilson ( 119751 ). In particular, configurations in the 
rapid and extreme rotation regimes exhibit a central toroidal 
structure, the volume of which increases with the value of the 
central rotation strength parameter. By making use of asymp- 
totic expansions of the density, mean velocity, and pressure 
tensor components for small radii, we found the condition 
for the existence of such central toroidal structure, as well as 
the condition for the maximum value of the central rotation 
strength parameter admitted by a configuration with a given 
concentration. 

- The differentially rotating models show a variety of realistic 
velocity dispersion profiles, characterized by the presence of 
pressure isotropy and radially-biased anisotropy in the cen- 
tral and intermediate regions, respectively. The kinematical 
behavior in the outer parts depends on the adopted trunca- 
tion prescription; in particular, the family which, in the non- 
rotating limit, reduces to the Wilson spheres is characterized 
by tangentially-biased anisotropy. This kinematical feature 
(rarely obtained in equilibrium models) is of great interest 
for two reasons: (i) Tangentially-biased pressure anisotropy 
is observed in the presence of internal rotation in globular 
clusters. For example, the full three-dimensional view of the 
velocity space of co Cen, obtained from proper motions and 
radial velocities measurements, has revealed that this ob- 
ject is characterized by significant rotation and tangential 
anisotropy in the outer parts (van de Ven et al. 20061). (ii) 
The dynamical evolution of a cluster in a tidal field is known 
to induce a rapid development of tangential anisotropy in the 
outer parts of the stellar system. In fact, if a cluster fills its 
Roche lobe and starts losing mass, there is a preferential loss 
of stars on radial orbits induced by the external tidal field 
at large radii, where tangential anisotropy in velocity space 
is thus established (Takahashi & Lee 2000; Baumgardt & 
Makino l2003l) . 

- The presence of differential rotation may induce nontriv- 
ial gradients in the line-of-sight velocity dispersion profile 
of a stellar system, even if the amount of rotation is mod- 
est. Therefore, this important physical ingredient should be 
taken into account properly. In this respect, dynamical stud- 
ies of globular clusters and other low-mass stellar systems 
by means of models based on the use of the Jeans equations 
are less satisfactory, because, at least in their most popular 
(nonrotating) application, they are used to reproduce varia- 
tions in the slope of the kinematical profile of a system only 
by means of a (sometimes significant) amount of pressure 
anisotropy. 

- As expected, differential rotation also induces nontrivial de- 
viations from spherical symmetry; in fact, the models are 
characterized by a great variety of (projected) ellipticity pro- 
files, dependent on the combined effect of concentration and 
central rotation strength. Configurations in the moderate ro- 
tation regime are characterized by realistic nonmonotonic 
ellipticity profiles (e.g., see Geyer et al. 119831) . while mod- 
els in the rapid rotation regime have monotonically decreas- 
ing profiles. To some extent, this morphological feature is 
complementary to that of the uniformly rotating models, in 
which the configurations are always characterized by mono- 
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tonically increasing ellipticity profiles. Interestingly, the be- 
havior of the ellipticity profiles is not necessarily correlated 
with the mean line-of-sight velocity profiles, that is configu- 
rations with a nonmonotonic mean velocity profile may have 
a monotonic ellipticity profile. In addition, the isophotes of 
the relevant surface density distribution tend to be character- 
ized by a "boxy" structure. 

- From a comparison of the equilibrium configurations result- 
ing from two options for the truncation prescription of the 
family of differentially rotating models, we confirm that the 
interplay between internal rotation, anisotropy in velocity 
space, and truncation in phase space is highly nontrivial. In 
fact, as also noted by Hunter J1977I) . the structure of the outer 
parts of a model is particularly sensitive, both from the mor- 
phological and the kinematical point of view, to the adopted 
truncation. One way to select the most appropriate truncation 
from the physical point of view will be to address the issue in 
the context of formation and evolution of the class of stellar 
systems under consideration (see also last item below). 

- Models in the moderate rotation regime seem to be partic- 
ularly appropriate for describing rotating globular clusters, 
since the relevant configurations are characterized by a num- 
ber of realistic properties, such as the presence of nonmono- 
tonic ellipticity profile, the behavior of surface density pro- 
file in the outer parts similar to the one associated with spher- 
ical Wilson models, the existence of pressure isotropy in 
the central regions and tangentially-biased anisotropy at the 
boundary, as well as realistic values of the ratio (vp)/crp. In 
a following paper, we plan to apply our family of differen- 
tially rotating models to selected Galactic globular clusters 
that show the presence of significant rotation, such as lo Cen 
and 47 Tuc. 

- Configurations with strong differential rotation, character- 
ized by the presence of a sizable central toroidal structure 
and by a off-centered peak of the surface brightness profiles, 
may be useful to shed light on the internal dynamics of the 
so-called "ring clusters". This class of object, originally ob- 
served in the Small Magellanic Cloud (Hill & Zaritsky 2006) 
and subsequently noted also in the Large Magellanic Cloud 
(Werchan & Zaritsky 1201 11 1, is characterized by a sizable 
dimple of the central surface brightness, resulting in an off- 
centered peaked density profile. A proper dynamical inter- 
pretation of these objects is currently missing. 

- The families of models illustrated in the present paper may 
also help to clarify the role of angular momentum in the for- 
mation and dynamical evolution of globular clusters. The re- 
sults of an extensive survey of N-body simulations, designed 
to study the dynamical stability and the long term evolution 
of the models described here, will be presented in follow-up 
papers (Varri et al. 12.01 It Varri et al. in preparation). 
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Appendix A: Spherical limit of rotating models 

The spherical nonrotating limit of the families of rotating mod- 
els considered in the present paper are defined by the following 
distribution functions 

f PT (E) = A&- aEa &- a(E - Eo) , (A.l) 
f K (E) = Ae- aE ° [e- a(£ - £o) - l] , (A.2) 

fwHE) = Ae- aEo [ e - fl(£ - £o) - 1 + a(E - Eq)] , (A.3) 

if E < Eq and fi(E) = otherwise (in the following, the in- 
dex i - 1,2,3 denotes the models defined by Eqs. ( 1A. U -( 1A.3I >. 
in the same order). In particular, /k(E) defines the King ( 1966) 
models and represents the spherical nonrotating limit of the fam- 
ily of rigidly rotating models defined by f' K {H) (see Eq. (Q}); 
fwr{E) and fpr(E) are the spherical isotropic limit of the Wilson 
d 19751 1 and Prendergast & Tomer (1 19701 1 models, and represent 
the spherical nonrotating limit of our differentially rotating mod- 
els defined by f* T (I) (see Eq. (ED) and fj T (I) (see Eq. d22]i), 
respectively. In the previous expressions, A, a, are positive con- 
stants, defining two dimensional scales, while Eg is the cut-off 
energy, which implies the existence of a truncation radius r tr 
for the spherical system. For all the families of models consid- 
ered here one important parameter is the central concentration, 
as measured by the depth of the dimensionless central potential 
well ¥ = tfr(Q) = a[Eo - O(0)] or by the parameter c = log(r /r ). 
Note that, if we compare three models (one for each family) 
having the same value of T, the corresponding values of c are 
not the same, that is the relevant dimensionless truncation radii 
are different (f, r increases from fpp(E), to /k(E), to fwr(E)). 
In other words, the structure of the outer parts of a spherical 
isotropic truncated model strictly depends on the truncation pre- 
scription; this property is particularly relevant for the interpreta- 
tion of the photometric profiles of globular clusters (see the sys- 
tematic comparison between spherical King and Wilson models 
performed by McLaughlin & van der Marel 2005). 

From the integration of the distribution functions in veloc- 
ity space, the corresponding intrinsic density distributions are 
recovered. Using the same dimensionless units introduced in 
the main text, we denote the dimensionless density profiles by 
Pi,s = Pi,s I A (for the definition of A, see Eq. (0), with 



here ip indicates the dimensionless escape energy, defined as in 
Eq.(ES. 

Similarly, the trace of the pressure tensor in dimensionless 
form (divided by a factor 3) can be written as = (a/A)p is , 
where 



(A.5) 



D,, P,, and £, are numerical coefficients resulting from the inte- 
gration in velocity space and are summarized in Tab. I A. 1 1 Note 
that the coefficient appearing as first argument of the incomplete 
gamma function in the pressure profile is related to the corre- 
sponding coefficient in the density profile, because they are, re- 
spectively, the second-order and zeroth-order moment of the dis- 
tribution function in velocity space. 

Table A.l. Coefficients for the spherical nonrotating models 



Model 


D 


P 


E 


Prendergast-Tomer 


3/2 


1 


3/2 


King 


1 


2/5 


5/2 


Wilson 


2/5 


4/35 


7/2 



Appendix B: Plain-truncated rotating models 

We summarize here the intrinsic properties of the family of mod- 
els defined by fp T (I) (see Eq. (l22l ). for a comparison with the 
corresponding quantities derived from the family of models de- 
fined by ft T (J) (see Eq. (TZTTi); a full description is given by Varri 
(120121) . The density profile in dimensionless form can be written 
as 



p PT (r, 0, i/0 = Pwrif, 6, I/O + i/r 



3/2 



~ 5 * 



5/2 . 



(B.l) 



the asymptotic behavior of the density profile in the outer parts 
(iff — > 0), with respect to the dimensionless escape energy, is 
given by: 

p PT (f, 6, M = if/' 2 + ^<A 5/2 (l + \xr- sin 2 flj + 0(if, 1/2 ) . (B.2) 

In the central region ( f — > 0), the expansion to second order in 
radius is 



p PT (f, 9, *P) = pPTfl + - 



9xr 2 sin 2 6»e' y r^, l pj + 



4 M2 dr 2 



r 2 + 0(?) 



(B.3) 



where the central value is ppr,o = 3/2e y(3/2,¥) = ppr^i^), 
consistent with the value of the corresponding nonrotating 
model. 

The mean velocity is in the azimuthal direction. In dimen- 
sionless units it is given by 



3e^ 



which, in the outer parts of the models reduces to 



dse~ s s J dttg(s,t,f,ff) , (B.4) 



p i ,s(i/r) = D i e*y(E i ,i/,) 



(A.4) ityprV, 0, = ^W l2 r sin 6 + 0(if, i/z ) . 



(B.5) 
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The expansion for f — > to first order in radius can be written as 



7(3/2, V ' 



(B.6) 



at variance with the family defined by f£ T (I), the dimensionless 
angular velocity depends also on the concentration parameter. 

As far as the pressure tensor is concerned, we recall that by 
construction p rr = pgg. The dimensionless radial component is 
given by 



,,7/2 



(B.7) 



which, expanded in 4> (as is appropriate for the outer parts), re- 
duces to 



(B.8) 



In the inner parts it can be approximated to second order in ra- 
dius, by the following expression 



PPT,n(f, 0, <¥) = ppTS) + i 



(B.9) 



where the central value is given by Ppt,o = e 4 'y(5/2, v P) = 
Ppt,sQ¥)i consistent with the value of the corresponding model 
in the limit of vanishing rotation. 

The azimuthal component of the pressure tensor can be writ- 
ten as 



pPT,w(r, e, (A) = Pw,m(?> 9 > M + 5 + -^f 11 + 

Pwt(v,/,)w T - Ppt{v4,) pt , 
which at the boundary is approximated by 



(B.10) 



7/2 I 



1 



■ X r l sm z 6) +0(if/ 9/2 ) 



35" \ 10 
Close to the center we find 



pPT44,{f, 0, TO = Ppts) + - \ 6% sin 6 



(B.ll) 
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2e 



4-7(5/2, TO 2 



7(3/2, T) 



3 » /3 

2 6r 2'* ^ 



r 2 +<9(r 4 ). (B.12) 



Therefore, the models in this family are characterized by 
pressure isotropy in the central region, radially-biased pressure 
anisotropy in the intermediate part, and pressure isotropy at the 
boundary (since, for both p rr and p^, the term of lower order in 
if/ is given by 2/5 t/^ 2 ), at variance with the family defined by 
f£ T (T) in which tangentially-biased anisotropy is present. 

Using the asymptotic expression of the density in the cen- 
tral regions recorded above, in this case, the condition for the 
existence of the central toroidal structure is given by 



X> 



7(1/2,40 
127(5/2, TO 



3 + 



C 2 /5 



2 \2 



1/2 



(B.13) 



where C2 is defined as in Eq. d38l >. By evaluating the sign of the 
velocity and pressure term in the radial component of the Jeans 
equation (see Eq. (l52l). we found that, in this case, the require- 
ment of positivity of the velocity term is just a necessary but not 
sufficient condition for the existence of the central toroidal struc- 
ture. In other words, in this family, configurations with angular 
velocity higher than the angular velocity associated with the cir- 
cular orbit of a single star but without a central toroidal structure 
can exist. The condition for maximally rotating configurations, 
which, in this case, is given by 



X< 



7d/2, T) /5 



1/2 



\C 2 \ 



(B.14) 



87(5/2, T) \2j 

completes the summary of the intrinsic properties of the family 
of plain-truncated differentially rotating models. 

Appendix C: The numerical procedure for the 
construction of differentially rotating models 

The construction of the equilibrium configurations for the fam- 
ilies of differentially rotating models, defined by Eqs. ( |2TT i and 
(l22l . has been performed by numerically solving the relevant 
Poisson equation as a nonlinear equation for the unknown po- 
tential ift. The code which implements the iteration procedure 
described in Sect. 13.21 starts with the calculation of the "seed 
solution", given by a selected spherical configuration from the 
family of models that represent the limit in the case of vanish- 
ing internal rotation of the family of interest (i.e., defined by 
Eqs. (IA.3b and ( lA.ll i. respectively). Such spherical solution is 
used, in the first step of the iteration, to evaluate (by means of 
a Double Gaussian Quadrature) the density distribution, given 
by pwt or ppj, on a spherical grid in the meridional plane, de- 
fined by (r,, 6j). The grid is linear in the radial coordinate and 
the angular positions are defined by 

„ x&l + 1) 



A L 



(CI) 



where j = 1, 2l max + 1. Typically, we used * 300 radial steps 
and we set l max = 21, in order to have sufficient accuracy to 
describe the complex morphologies of the configurations in the 
rapid and extreme rotation regimes. The discrete direct and in- 
verse Legendre transforms, required at each step of the iteration 
for the calculation of the density and potential coefficients, de- 
fined by Eqs. (t26b-(f25b. are performed (up to the order l max ) by 
means of a package based on S2kit 1 .0 by Kostelec & Rockmore 

(2004) , which makes use of FFTW 3.2.1 by Frigo & Johnston 

(2005) . The Cauchy problems for the potential radial functions 
expressed in Eqs. d3T1l-(l32b are evaluated by means of a numeri- 
cal integration with Romberg's rule. The convergence condition 
for the solution at the n-th step of the iteration is formally defined 
as 



11 < 10- 3 



(C2) 



for every where if/r9 = ^"\ri,6j); about 10 (25) iteration 
steps are needed for the construction of configurations charac- 
terized by low (high) values of %■ The accuracy of the solutions 
found with our code has been checked by the following tests: 
(i) the virial theorem is satisfied with accuracy of the order of 
10~ 4 or better; (2) the radial component of the Jeans equation is 
satisfied with the accuracy of the order of 10~ 3 or better; (3) the 
asymptotic behaviors, both in the central and in the outer parts, 
of all the moments of the distribution function are confirmed. 
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